/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */


#include "bsr_common.h"
#include "parbsr_laplacian27pt_generator.h"

namespace YHAMG
{

ParBSR_Laplacian27pt_Generator::ParBSR_Laplacian27pt_Generator()
	: nx(50),
	ny(50),
	nz(50),
	bx(2),
	by(2),
	bz(2),
	Px(1),
	Py(1),
	Pz(1)
{
}

ParBSR_Laplacian27pt_Generator::ParBSR_Laplacian27pt_Generator(MPI_Comm comm)
	: ParBSR_Generator(comm),
	nx(50),
	ny(50),
	nz(50),
	bx(2),
	by(2),
	bz(2),
	Px(1),
	Py(1),
	Pz(1)
{
}

ParBSR_Laplacian27pt_Generator::ParBSR_Laplacian27pt_Generator(MPI_Comm comm, int _nx, int _ny, int _nz, int _bx, int _by, int _bz, int _Px, int _Py, int _Pz)
	: ParBSR_Generator(comm),
	nx(_nx),
	ny(_ny),
	nz(_nz),
	bx(_bx),
	by(_by),
	bz(_bz),
	Px(_Px),
	Py(_Py),
	Pz(_Pz)
{
}

#define _N(ix, iy, iz) ((ix) * ny * nz + (iy) * nz + (iz))
#define _B(jx, jy, jz) ((jx) * by * bz + (jy) * bz + (jz))
#define _P(rx, ry, rz) ((rx) * Py * Pz + (ry) * Pz + (rz))

void ParBSR_Laplacian27pt_Generator::operator()(ParBSR_Matrix& A) const
{
	int rank;
	MPI_Comm_rank(comm, &rank);

	double c = 1.0;
	double d = -1.0;
	if (nx * Px * bx > 1) d *= 3;
	if (ny * Py * bx > 1) d *= 3;
	if (nz * Pz * bx > 1) d *= 3;
	d += 1.0;

	int size = nx * ny * nz;
	int block_size = bx * by * bz;
	int block_nnz = block_size * block_size;

	int* LocLoc_rowptr = new int[size + 1];

	LocLoc_rowptr[0] = 0;
	for (int ix = 0; ix < nx; ++ix)
	{
		for (int iy = 0; iy < ny; ++iy)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				int cnt = 1;
				if (ix > 0)
				{
					++cnt;
					if (iy > 0)
					{
						++cnt;
						if (iz > 0) ++cnt;
						if (iz < nz - 1) ++cnt;
					}
					if (iy < ny - 1)
					{
						++cnt;
						if (iz > 0) ++cnt;
						if (iz < nz - 1) ++cnt;
					}
					if (iz > 0) ++cnt;
					if (iz < nz - 1) ++cnt;
				}

				if (ix < nx - 1)
				{
					++cnt;
					if (iy > 0)
					{
						++cnt;
						if (iz > 0) ++cnt;
						if (iz < nz - 1) ++cnt;
					}
					if (iy < ny - 1)
					{
						++cnt;
						if (iz > 0) ++cnt;
						if (iz < nz - 1) ++cnt;
					}
					if (iz > 0) ++cnt;
					if (iz < nz - 1) ++cnt;
				}

				if (iy > 0)
				{
					++cnt;
					if (iz > 0) ++cnt;
					if (iz < nz - 1) ++cnt;
				}

				if (iy < ny - 1)
				{
					++cnt;
					if (iz > 0) ++cnt;
					if (iz < nz - 1) ++cnt;
				}

				if (iz > 0) ++cnt;
				if (iz < nz - 1) ++cnt;

				LocLoc_rowptr[_N(ix, iy, iz) + 1] = cnt;
			}
		}
	}

	for (int i = 0; i < size; ++i)
		LocLoc_rowptr[i + 1] += LocLoc_rowptr[i];

	int* LocLoc_colind = new int[LocLoc_rowptr[size]];
	double* LocLoc_values = new double[LocLoc_rowptr[size] * block_nnz];

	for (int ix = 0; ix < nx; ++ix)
	{
		for (int iy = 0; iy < ny; ++iy)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				int j = LocLoc_rowptr[_N(ix, iy, iz)];

				if (ix > 0)
				{
					if (iy > 0)
					{
						if (iz > 0)
						{
							LocLoc_colind[j] = _N(ix - 1, iy - 1, iz - 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
							++j;
						}

						LocLoc_colind[j] = _N(ix - 1, iy - 1, iz);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jz = 0; jz < bz; ++jz)
						{
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(bx - 1, by - 1, jz - 1) * block_size + _B(0, 0, jz)] = c;
							LocLoc_values[j * block_nnz + _B(bx - 1, by - 1, jz) * block_size + _B(0, 0, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(bx - 1, by - 1, jz + 1) * block_size + _B(0, 0, jz)] = c;
						}

						++j;

						if (iz < nz - 1)
						{
							LocLoc_colind[j] = _N(ix - 1, iy - 1, iz + 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
							++j;
						}
					}

					if (iz > 0)
					{
						LocLoc_colind[j] = _N(ix - 1, iy, iz - 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
			
						for (int jy = 0; jy < by; ++jy)
						{
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy - 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
							LocLoc_values[j * block_nnz + _B(bx - 1, jy, bz - 1) * block_size + _B(0, jy, 0)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy + 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
						}

						++j;
					}

					LocLoc_colind[j] = _N(ix - 1, iy, iz);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jy = 0; jy < by; ++jy)
					{
						for (int jz = 0; jz < bz; ++jz)
						{
							if (jy > 0)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(bx - 1, jy - 1, jz - 1) * block_size + _B(0, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(bx - 1, jy - 1, jz) * block_size + _B(0, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(bx - 1, jy - 1, jz + 1) * block_size + _B(0, jy, jz)] = c;
							}
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy, jz - 1) * block_size + _B(0, jy, jz)] = c;
							LocLoc_values[j * block_nnz + _B(bx - 1, jy, jz) * block_size + _B(0, jy, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy, jz + 1) * block_size + _B(0, jy, jz)] = c;
							if (jy < by - 1)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(bx - 1, jy + 1, jz - 1) * block_size + _B(0, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(bx - 1, jy + 1, jz) * block_size + _B(0, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(bx - 1, jy + 1, jz + 1) * block_size + _B(0, jy, jz)] = c;
							}
						}
					}

					++j;

					if (iz < nz - 1)
					{
						LocLoc_colind[j] = _N(ix - 1, iy, iz + 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
			
						for (int jy = 0; jy < by; ++jy)
						{
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy - 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
							LocLoc_values[j * block_nnz + _B(bx - 1, jy, 0) * block_size + _B(0, jy, bz - 1)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(bx - 1, jy + 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
						}

						++j;
					}

					if (iy < ny - 1)
					{
						if (iz > 0)
						{
							LocLoc_colind[j] = _N(ix - 1, iy + 1, iz - 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
							++j;
						}

						LocLoc_colind[j] = _N(ix - 1, iy + 1, iz);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jz = 0; jz < bz; ++jz)
						{
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(bx - 1, 0, jz - 1) * block_size + _B(0, by - 1, jz)] = c;
							LocLoc_values[j * block_nnz + _B(bx - 1, 0, jz) * block_size + _B(0, by - 1, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(bx - 1, 0, jz + 1) * block_size + _B(0, by - 1, jz)] = c;
						}

						++j;

						if (iz < nz - 1)
						{
							LocLoc_colind[j] = _N(ix - 1, iy + 1, iz + 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
							++j;
						}
					}
				}

				if (iy > 0)
				{
					if (iz > 0)
					{
						LocLoc_colind[j] = _N(ix, iy - 1, iz - 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jx = 0; jx < bx; ++jx)
						{
							if (jx > 0)
								LocLoc_values[j * block_nnz + _B(jx - 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
							LocLoc_values[j * block_nnz + _B(jx, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
							if (jx < bx - 1)
								LocLoc_values[j * block_nnz + _B(jx + 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
						}

						++j;
					}

					LocLoc_colind[j] = _N(ix, iy - 1, iz);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jx = 0; jx < bx; ++jx)
					{
						for (int jz = 0; jz < bz; ++jz)
						{
							if (jx > 0)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx - 1, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx - 1, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx - 1, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
							}
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(jx, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
							LocLoc_values[j * block_nnz + _B(jx, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(jx, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
							if (jx < bx - 1)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx + 1, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx + 1, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx + 1, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
							}
						}
					}

					++j;

					if (iz < nz - 1)
					{
						LocLoc_colind[j] = _N(ix, iy - 1, iz + 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jx = 0; jx < bx; ++jx)
						{
							if (jx > 0)
								LocLoc_values[j * block_nnz + _B(jx - 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
							LocLoc_values[j * block_nnz + _B(jx, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
							if (jx < bx - 1)
								LocLoc_values[j * block_nnz + _B(jx + 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
						}

						++j;
					}
				}

				if (iz > 0)
				{
					LocLoc_colind[j] = _N(ix, iy, iz - 1);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jx = 0; jx < bx; ++jx)
					{
						for (int jy = 0; jy < by; ++jy)
						{
							if (jx > 0)
							{
								if (jy > 0)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
								LocLoc_values[j * block_nnz + _B(jx - 1, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
								if (jy < by - 1)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							}
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(jx, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							LocLoc_values[j * block_nnz + _B(jx, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(jx, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							if (jx < bx - 1)
							{
								if (jy > 0)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
								LocLoc_values[j * block_nnz + _B(jx + 1, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
								if (jy < by - 1)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							}
						}
					}

					++j;
				}

				LocLoc_colind[j] = _N(ix, iy, iz);
				BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

				for (int jx = 0; jx < bx; ++jx)
				{
					for (int jy = 0; jy < by; ++jy)
					{
						for (int jz = 0; jz < bz; ++jz)
						{
							if (jx > 0)
							{
								if (jy > 0)
								{
									if (jz > 0)
										LocLoc_values[j * block_nnz + _B(jx - 1, jy - 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
									LocLoc_values[j * block_nnz + _B(jx - 1, jy - 1, jz) * block_size + _B(jx, jy, jz)] = c;
									if (jz < bz - 1)
										LocLoc_values[j * block_nnz + _B(jx - 1, jy - 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								}
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy, jz - 1) * block_size + _B(jx, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx - 1, jy, jz) * block_size + _B(jx, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								if (jy < by - 1)
								{
									if (jz > 0)
										LocLoc_values[j * block_nnz + _B(jx - 1, jy + 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
									LocLoc_values[j * block_nnz + _B(jx - 1, jy + 1, jz) * block_size + _B(jx, jy, jz)] = c;
									if (jz < bz - 1)
										LocLoc_values[j * block_nnz + _B(jx - 1, jy + 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								}
							}

							if (jy > 0)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx, jy - 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx, jy - 1, jz) * block_size + _B(jx, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx, jy - 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
							}
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(jx, jy, jz - 1) * block_size + _B(jx, jy, jz)] = c;
							LocLoc_values[j * block_nnz + _B(jx, jy, jz) * block_size + _B(jx, jy, jz)] = d;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(jx, jy, jz + 1) * block_size + _B(jx, jy, jz)] = c;
							if (jy < by - 1)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx, jy + 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx, jy + 1, jz) * block_size + _B(jx, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx, jy + 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
							}
							
							if (jx < bx - 1)
							{
								if (jy > 0)
								{
									if (jz > 0)
										LocLoc_values[j * block_nnz + _B(jx + 1, jy - 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
									LocLoc_values[j * block_nnz + _B(jx + 1, jy - 1, jz) * block_size + _B(jx, jy, jz)] = c;
									if (jz < bz - 1)
										LocLoc_values[j * block_nnz + _B(jx + 1, jy - 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								}
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy, jz - 1) * block_size + _B(jx, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx + 1, jy, jz) * block_size + _B(jx, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								if (jy < by - 1)
								{
									if (jz > 0)
										LocLoc_values[j * block_nnz + _B(jx + 1, jy + 1, jz - 1) * block_size + _B(jx, jy, jz)] = c;
									LocLoc_values[j * block_nnz + _B(jx + 1, jy + 1, jz) * block_size + _B(jx, jy, jz)] = c;
									if (jz < bz - 1)
										LocLoc_values[j * block_nnz + _B(jx + 1, jy + 1, jz + 1) * block_size + _B(jx, jy, jz)] = c;
								}
							}
						}
					}
				}

				++j;

				if (iz < nz - 1)
				{
					LocLoc_colind[j] = _N(ix, iy, iz + 1);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jx = 0; jx < bx; ++jx)
					{
						for (int jy = 0; jy < by; ++jy)
						{
							if (jx > 0)
							{
								if (jy > 0)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
								LocLoc_values[j * block_nnz + _B(jx - 1, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
								if (jy < by - 1)
									LocLoc_values[j * block_nnz + _B(jx - 1, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							}
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(jx, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							LocLoc_values[j * block_nnz + _B(jx, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(jx, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							if (jx < bx - 1)
							{
								if (jy > 0)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
								LocLoc_values[j * block_nnz + _B(jx + 1, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
								if (jy < by - 1)
									LocLoc_values[j * block_nnz + _B(jx + 1, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							}
						}
					}

					++j;
				}

				if (iy < ny - 1)
				{
					if (iz > 0)
					{
						LocLoc_colind[j] = _N(ix, iy + 1, iz - 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jx = 0; jx < bx; ++jx)
						{
							if (jx > 0)
								LocLoc_values[j * block_nnz + _B(jx - 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
							LocLoc_values[j * block_nnz + _B(jx, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
							if (jx < bx - 1)
								LocLoc_values[j * block_nnz + _B(jx + 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
						}

						++j;
					}

					LocLoc_colind[j] = _N(ix, iy + 1, iz);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jx = 0; jx < bx; ++jx)
					{
						for (int jz = 0; jz < bz; ++jz)
						{
							if (jx > 0)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx - 1, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx - 1, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx - 1, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
							}
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(jx, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
							LocLoc_values[j * block_nnz + _B(jx, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(jx, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
							if (jx < bx - 1)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(jx + 1, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
								LocLoc_values[j * block_nnz + _B(jx + 1, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(jx + 1, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
							}
						}
					}

					++j;

					if (iz < nz - 1)
					{
						LocLoc_colind[j] = _N(ix, iy + 1, iz + 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jx = 0; jx < bx; ++jx)
						{
							if (jx > 0)
								LocLoc_values[j * block_nnz + _B(jx - 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
							LocLoc_values[j * block_nnz + _B(jx, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
							if (jx < bx - 1)
								LocLoc_values[j * block_nnz + _B(jx + 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
						}

						++j;
					}
				}

				if (ix < nx - 1)
				{
					if (iy > 0)
					{
						if (iz > 0)
						{
							LocLoc_colind[j] = _N(ix + 1, iy - 1, iz - 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
							++j;
						}

						LocLoc_colind[j] = _N(ix + 1, iy - 1, iz);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jz = 0; jz < bz; ++jz)
						{
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(0, by - 1, jz - 1) * block_size + _B(bx - 1, 0, jz)] = c;
							LocLoc_values[j * block_nnz + _B(0, by - 1, jz) * block_size + _B(bx - 1, 0, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(0, by - 1, jz + 1) * block_size + _B(bx - 1, 0, jz)] = c;
						}

						++j;

						if (iz < nz - 1)
						{
							LocLoc_colind[j] = _N(ix + 1, iy - 1, iz + 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
							++j;
						}
					}

					if (iz > 0)
					{
						LocLoc_colind[j] = _N(ix + 1, iy, iz - 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
			
						for (int jy = 0; jy < by; ++jy)
						{
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(0, jy - 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
							LocLoc_values[j * block_nnz + _B(0, jy, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(0, jy + 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
						}

						++j;
					}

					LocLoc_colind[j] = _N(ix + 1, iy, iz);
					BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

					for (int jy = 0; jy < by; ++jy)
					{
						for (int jz = 0; jz < bz; ++jz)
						{
							if (jy > 0)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(0, jy - 1, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(0, jy - 1, jz) * block_size + _B(bx - 1, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(0, jy - 1, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
							}
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(0, jy, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
							LocLoc_values[j * block_nnz + _B(0, jy, jz) * block_size + _B(bx - 1, jy, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(0, jy, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
							if (jy < by - 1)
							{
								if (jz > 0)
									LocLoc_values[j * block_nnz + _B(0, jy + 1, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
								LocLoc_values[j * block_nnz + _B(0, jy + 1, jz) * block_size + _B(bx - 1, jy, jz)] = c;
								if (jz < bz - 1)
									LocLoc_values[j * block_nnz + _B(0, jy + 1, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
							}
						}
					}

					++j;

					if (iz < nz - 1)
					{
						LocLoc_colind[j] = _N(ix + 1, iy, iz + 1);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
			
						for (int jy = 0; jy < by; ++jy)
						{
							if (jy > 0)
								LocLoc_values[j * block_nnz + _B(0, jy - 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
							LocLoc_values[j * block_nnz + _B(0, jy, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
							if (jy < by - 1)
								LocLoc_values[j * block_nnz + _B(0, jy + 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
						}

						++j;
					}

					if (iy < ny - 1)
					{
						if (iz > 0)
						{
							LocLoc_colind[j] = _N(ix + 1, iy + 1, iz - 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
							++j;
						}

						LocLoc_colind[j] = _N(ix + 1, iy + 1, iz);
						BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);

						for (int jz = 0; jz < bz; ++jz)
						{
							if (jz > 0)
								LocLoc_values[j * block_nnz + _B(0, 0, jz - 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
							LocLoc_values[j * block_nnz + _B(0, 0, jz) * block_size + _B(bx - 1, by - 1, jz)] = c;
							if (jz < bz - 1)
								LocLoc_values[j * block_nnz + _B(0, 0, jz + 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
						}

						++j;

						if (iz < nz - 1)
						{
							LocLoc_colind[j] = _N(ix + 1, iy + 1, iz + 1);
							BSRBlockFill(block_size, 0.0, LocLoc_values + j * block_nnz);
							LocLoc_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
							++j;
						}
					}
				}
			}
		}
	}

	int px = rank / (Py * Pz);
	int py = rank / Pz - px * Py;
	int pz = rank - px * Py * Pz - py * Pz;

	int* LocExt_rowptr = new int[size + 1];
	for (int i = 0; i <= size; ++i)
		LocExt_rowptr[i] = 0;

	int num_neighbors = 0;
	int* neighbor_ranks = new int[26];
	int* Ext_ptr = new int[27];

	Ext_ptr[0] = 0;

	if (px > 0)
	{
		if (py > 0)
		{
			if (pz > 0)
			{
				++LocExt_rowptr[_N(0, 0, 0)];
				neighbor_ranks[num_neighbors++] = _P(px - 1, py - 1, pz - 1);
				Ext_ptr[num_neighbors] = 1;
			}

			for (int iz = 0; iz < nz; ++iz)
			{
				if (iz > 0) ++LocExt_rowptr[_N(0, 0, iz)];
				++LocExt_rowptr[_N(0, 0, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(0, 0, iz)];
			}
			neighbor_ranks[num_neighbors++] = _P(px - 1, py - 1, pz);
			Ext_ptr[num_neighbors] = nz;

			if (pz < Pz - 1)
			{
				++LocExt_rowptr[_N(0, 0, nz - 1)];
				neighbor_ranks[num_neighbors++] = _P(px - 1, py - 1, pz + 1);
				Ext_ptr[num_neighbors] = 1;
			}
		}

		if (pz > 0)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (iy > 0) ++LocExt_rowptr[_N(0, iy, 0)];
				++LocExt_rowptr[_N(0, iy, 0)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(0, iy, 0)];
			}
			neighbor_ranks[num_neighbors++] = _P(px - 1, py, pz - 1);
			Ext_ptr[num_neighbors] = ny;
		}

		for (int iy = 0; iy < ny; ++iy)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				if (iy > 0)
				{
					if (iz > 0) ++LocExt_rowptr[_N(0, iy, iz)];
					++LocExt_rowptr[_N(0, iy, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(0, iy, iz)];
				}
				if (iz > 0) ++LocExt_rowptr[_N(0, iy, iz)];
				++LocExt_rowptr[_N(0, iy, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(0, iy, iz)];
				if (iy < ny - 1)
				{
					if (iz > 0) ++LocExt_rowptr[_N(0, iy, iz)];
					++LocExt_rowptr[_N(0, iy, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(0, iy, iz)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px - 1, py, pz);
		Ext_ptr[num_neighbors] = ny * nz;

		if (pz < Pz - 1)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (iy > 0) ++LocExt_rowptr[_N(0, iy, nz - 1)];
				++LocExt_rowptr[_N(0, iy, nz - 1)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(0, iy, nz - 1)];
			}
			neighbor_ranks[num_neighbors++] = _P(px - 1, py, pz + 1);
			Ext_ptr[num_neighbors] = ny;
		}

		if (py < Py - 1)
		{
			if (pz > 0)
			{
				++LocExt_rowptr[_N(0, ny - 1, 0)];
				neighbor_ranks[num_neighbors++] = _P(px - 1, py + 1, pz - 1);
				Ext_ptr[num_neighbors] = 1;
			}

			for (int iz = 0; iz < nz; ++iz)
			{
				if (iz > 0) ++LocExt_rowptr[_N(0, ny - 1, iz)];
				++LocExt_rowptr[_N(0, ny - 1, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(0, ny - 1, iz)];
			}
			neighbor_ranks[num_neighbors++] = _P(px - 1, py + 1, pz);
			Ext_ptr[num_neighbors] = nz;

			if (pz < Pz - 1)
			{
				++LocExt_rowptr[_N(0, ny - 1, nz - 1)];
				neighbor_ranks[num_neighbors++] = _P(px - 1, py + 1, pz + 1);
				Ext_ptr[num_neighbors] = 1;
			}
		}
	}

	if (py > 0)
	{
		if (pz > 0)
		{
			for (int ix = 0; ix < nx; ++ix)
			{
				if (ix > 0) ++LocExt_rowptr[_N(ix, 0, 0)];
				++LocExt_rowptr[_N(ix, 0, 0)];
				if (ix < nx - 1) ++LocExt_rowptr[_N(ix, 0, 0)];
			}
			neighbor_ranks[num_neighbors++] = _P(px, py - 1, pz - 1);
			Ext_ptr[num_neighbors] = nx;
		}

		for (int ix = 0; ix < nx; ++ix)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				if (ix > 0)
				{
					if (iz > 0) ++LocExt_rowptr[_N(ix, 0, iz)];
					++LocExt_rowptr[_N(ix, 0, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(ix, 0, iz)];
				}
				if (iz > 0) ++LocExt_rowptr[_N(ix, 0, iz)];
				++LocExt_rowptr[_N(ix, 0, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(ix, 0, iz)];
				if (ix < nx - 1)
				{
					if (iz > 0) ++LocExt_rowptr[_N(ix, 0, iz)];
					++LocExt_rowptr[_N(ix, 0, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(ix, 0, iz)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px, py - 1, pz);
		Ext_ptr[num_neighbors] = nx * nz;

		if (pz < Pz - 1)
		{
			for (int ix = 0; ix < nx; ++ix)
			{
				if (ix > 0) ++LocExt_rowptr[_N(ix, 0, nz - 1)];
				++LocExt_rowptr[_N(ix, 0, nz - 1)];
				if (ix < nx - 1) ++LocExt_rowptr[_N(ix, 0, nz - 1)];
			}
			neighbor_ranks[num_neighbors++] = _P(px, py - 1, pz + 1);
			Ext_ptr[num_neighbors] = nx;
		}
	}

	if (pz > 0)
	{
		for (int ix = 0; ix < nx; ++ix)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (ix > 0)
				{
					if (iy > 0) ++LocExt_rowptr[_N(ix, iy, 0)];
					++LocExt_rowptr[_N(ix, iy, 0)];
					if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, 0)];
				}
				if (iy > 0) ++LocExt_rowptr[_N(ix, iy, 0)];
				++LocExt_rowptr[_N(ix, iy, 0)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, 0)];
				if (ix < nx - 1)
				{
					if (iy > 0) ++LocExt_rowptr[_N(ix, iy, 0)];
					++LocExt_rowptr[_N(ix, iy, 0)];
					if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, 0)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px, py, pz - 1);
		Ext_ptr[num_neighbors] = nx * ny;
	}

	if (pz < Pz - 1)
	{
		for (int ix = 0; ix < nx; ++ix)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (ix > 0)
				{
					if (iy > 0) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
					++LocExt_rowptr[_N(ix, iy, nz - 1)];
					if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
				}
				if (iy > 0) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
				++LocExt_rowptr[_N(ix, iy, nz - 1)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
				if (ix < nx - 1)
				{
					if (iy > 0) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
					++LocExt_rowptr[_N(ix, iy, nz - 1)];
					if (iy < ny - 1) ++LocExt_rowptr[_N(ix, iy, nz - 1)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px, py, pz + 1);
		Ext_ptr[num_neighbors] = nx * ny;
	}

	if (py < Py - 1)
	{
		if (pz > 0)
		{
			for (int ix = 0; ix < nx; ++ix)
			{
				if (ix > 0) ++LocExt_rowptr[_N(ix, ny - 1, 0)];
				++LocExt_rowptr[_N(ix, ny - 1, 0)];
				if (ix < nx - 1) ++LocExt_rowptr[_N(ix, ny - 1, 0)];
			}
			neighbor_ranks[num_neighbors++] = _P(px, py + 1, pz - 1);
			Ext_ptr[num_neighbors] = nx;
		}

		for (int ix = 0; ix < nx; ++ix)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				if (ix > 0)
				{
					if (iz > 0) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
					++LocExt_rowptr[_N(ix, ny - 1, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
				}
				if (iz > 0) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
				++LocExt_rowptr[_N(ix, ny - 1, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
				if (ix < nx - 1)
				{
					if (iz > 0) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
					++LocExt_rowptr[_N(ix, ny - 1, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(ix, ny - 1, iz)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px, py + 1, pz);
		Ext_ptr[num_neighbors] = nx * nz;

		if (pz < Pz - 1)
		{
			for (int ix = 0; ix < nx; ++ix)
			{
				if (ix > 0) ++LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
				++LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
				if (ix < nx - 1) ++LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
			}
			neighbor_ranks[num_neighbors++] = _P(px, py + 1, pz + 1);
			Ext_ptr[num_neighbors] = nx;
		}
	}

	if (px < Px - 1)
	{
		if (py > 0)
		{
			if (pz > 0)
			{
				++LocExt_rowptr[_N(nx - 1, 0, 0)];
				neighbor_ranks[num_neighbors++] = _P(px + 1, py - 1, pz - 1);
				Ext_ptr[num_neighbors] = 1;
			}

			for (int iz = 0; iz < nz; ++iz)
			{
				if (iz > 0) ++LocExt_rowptr[_N(nx - 1, 0, iz)];
				++LocExt_rowptr[_N(nx - 1, 0, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(nx - 1, 0, iz)];
			}
			neighbor_ranks[num_neighbors++] = _P(px + 1, py - 1, pz);
			Ext_ptr[num_neighbors] = nz;

			if (pz < Pz - 1)
			{
				++LocExt_rowptr[_N(nx - 1, 0, nz - 1) ];
				neighbor_ranks[num_neighbors++] = _P(px + 1, py - 1, pz + 1);
				Ext_ptr[num_neighbors] = 1;
			}
		}

		if (pz > 0)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (iy > 0) ++LocExt_rowptr[_N(nx - 1, iy, 0)];
				++LocExt_rowptr[_N(nx - 1, iy, 0)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(nx - 1, iy, 0)];
			}
			neighbor_ranks[num_neighbors++] = _P(px + 1, py, pz - 1);
			Ext_ptr[num_neighbors] = ny;
		}

		for (int iy = 0; iy < ny; ++iy)
		{
			for (int iz = 0; iz < nz; ++iz)
			{
				if (iy > 0)
				{
					if (iz > 0) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
					++LocExt_rowptr[_N(nx - 1, iy, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
				}
				if (iz > 0) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
				++LocExt_rowptr[_N(nx - 1, iy, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
				if (iy < ny - 1)
				{
					if (iz > 0) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
					++LocExt_rowptr[_N(nx - 1, iy, iz)];
					if (iz < nz - 1) ++LocExt_rowptr[_N(nx - 1, iy, iz)];
				}
			}
		}
		neighbor_ranks[num_neighbors++] = _P(px + 1, py, pz);
		Ext_ptr[num_neighbors] = ny * nz;

		if (pz < Pz - 1)
		{
			for (int iy = 0; iy < ny; ++iy)
			{
				if (iy > 0) ++LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
				++LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
				if (iy < ny - 1) ++LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
			}
			neighbor_ranks[num_neighbors++] = _P(px + 1, py, pz + 1);
			Ext_ptr[num_neighbors] = ny;
		}

		if (py < Py - 1)
		{
			if (pz > 0)
			{
				++LocExt_rowptr[_N(nx - 1, ny - 1, 0)];
				neighbor_ranks[num_neighbors++] = _P(px + 1, py + 1, pz - 1);
				Ext_ptr[num_neighbors] = 1;
			}

			for (int iz = 0; iz < nz; ++iz)
			{
				if (iz > 0) ++LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
				++LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
				if (iz < nz - 1) ++LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
			}
			neighbor_ranks[num_neighbors++] = _P(px + 1, py + 1, pz);
			Ext_ptr[num_neighbors] = nz;

			if (pz < Pz - 1)
			{
				++LocExt_rowptr[_N(nx - 1, ny - 1, nz - 1)];
				neighbor_ranks[num_neighbors++] = _P(px + 1, py + 1, pz + 1);
				Ext_ptr[num_neighbors] = 1;
			}
		}
	}

	for (int i = 0; i < size; ++i)
		LocExt_rowptr[i + 1] += LocExt_rowptr[i];

	int* LocExt_colind = new int[LocExt_rowptr[size]];
	double* LocExt_values = new double[LocExt_rowptr[size] * block_nnz];

	for (int r = 0; r < num_neighbors; ++r)
		Ext_ptr[r + 1] += Ext_ptr[r];

	int* Ext_map = new int[Ext_ptr[num_neighbors]];

	int neighbor_index = num_neighbors;

	if (px < Px - 1)
	{
		if (py < Py - 1)
		{
			if (pz < Pz - 1)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(0, 0, 0);
				int j = --LocExt_rowptr[_N(nx - 1, ny - 1, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
			}

			int t = Ext_ptr[neighbor_index--];
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(0, 0, iz);
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jz = 0; jz < bz; ++jz)
				{
					if (jz > 0)
						LocExt_values[j * block_nnz + _B(0, 0, jz - 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
					LocExt_values[j * block_nnz + _B(0, 0, jz) * block_size + _B(bx - 1, by - 1, jz)] = c;
					if (jz < bz - 1)
						LocExt_values[j * block_nnz + _B(0, 0, jz + 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(nx - 1, ny - 1, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
				}
			}

			if (pz > 0)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(0, 0, nz - 1);
				int j = --LocExt_rowptr[_N(nx - 1, ny - 1, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
			}
		}

		if (pz < Pz - 1)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(0, iy, 0);
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					if (jy > 0)
						LocExt_values[j * block_nnz + _B(0, jy - 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
					LocExt_values[j * block_nnz + _B(0, jy, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
					if (jy < by - 1)
						LocExt_values[j * block_nnz + _B(0, jy + 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, nz - 1)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
				}
			}
		}

		int t = Ext_ptr[neighbor_index--];
		for (int iy = ny - 1; iy >= 0; --iy)
		{
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(0, iy, iz);
				if (iy < ny - 1)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
						LocExt_colind[j] = t + nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
					LocExt_colind[j] = t + nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(0, 0, jz - 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
						LocExt_values[j * block_nnz + _B(0, 0, jz) * block_size + _B(bx - 1, by - 1, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(0, 0, jz + 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
						LocExt_colind[j] = t + nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
					}
				}
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(0, jy - 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(0, jy, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(0, jy + 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jy > 0)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(0, jy - 1, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
							LocExt_values[j * block_nnz + _B(0, jy - 1, jz) * block_size + _B(bx - 1, jy, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(0, jy - 1, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
						}
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(0, jy, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
						LocExt_values[j * block_nnz + _B(0, jy, jz) * block_size + _B(bx - 1, jy, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(0, jy, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
						if (jy < by - 1)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(0, jy + 1, jz - 1) * block_size + _B(bx - 1, jy, jz)] = c;
							LocExt_values[j * block_nnz + _B(0, jy + 1, jz) * block_size + _B(bx - 1, jy, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(0, jy + 1, jz + 1) * block_size + _B(bx - 1, jy, jz)] = c;
						}
					}
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(0, jy - 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
						LocExt_values[j * block_nnz + _B(0, jy, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(0, jy + 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
					}
				}
				if (iy > 0)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
						LocExt_colind[j] = t - nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
					LocExt_colind[j] = t - nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(0, by - 1, jz - 1) * block_size + _B(bx - 1, 0, jz)] = c;
						LocExt_values[j * block_nnz + _B(0, by - 1, jz) * block_size + _B(bx - 1, 0, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(0, by - 1, jz + 1) * block_size + _B(bx - 1, 0, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(nx - 1, iy, iz)];
						LocExt_colind[j] = t - nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
					}
				}
			}
		}

		if (pz > 0)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(0, iy, nz - 1);
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, 0)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
				}
				int j = --LocExt_rowptr[_N(nx - 1, iy, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					if (jy > 0)
						LocExt_values[j * block_nnz + _B(0, jy - 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
					LocExt_values[j * block_nnz + _B(0, jy, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
					if (jy < by - 1)
						LocExt_values[j * block_nnz + _B(0, jy + 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(nx - 1, iy, 0)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
				}
			}
		}

		if (py > 0)
		{
			if (pz < Pz - 1)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(0, ny - 1, 0);
				int j = --LocExt_rowptr[_N(nx - 1, 0, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
			}

			int t = Ext_ptr[neighbor_index--];
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(0, ny - 1, iz);
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(nx - 1, 0, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(nx - 1, 0, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jz = 0; jz < bz; ++jz)
				{
					if (jz > 0)
						LocExt_values[j * block_nnz + _B(0, by - 1, jz - 1) * block_size + _B(bx - 1, 0, jz)] = c;
					LocExt_values[j * block_nnz + _B(0, by - 1, jz) * block_size + _B(bx - 1, 0, jz)] = c;
					if (jz < bz - 1)
						LocExt_values[j * block_nnz + _B(0, by - 1, jz + 1) * block_size + _B(bx - 1, 0, jz)] = c;
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(nx - 1, 0, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
				}
			}

			if (pz > 0)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(0, ny - 1, nz - 1);
				int j = --LocExt_rowptr[_N(nx - 1, 0, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
			}
		}
	}

	if (py < Py - 1)
	{
		if (pz < Pz - 1)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int ix = nx - 1; ix >= 0; --ix)
			{
				Ext_map[--t] = _N(ix, 0, 0);
				if (ix < nx - 1)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					if (jx > 0)
						LocExt_values[j * block_nnz + _B(jx - 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(jx, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
					if (jx < bx - 1)
						LocExt_values[j * block_nnz + _B(jx + 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
				}
				if (ix > 0)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, nz - 1)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
				}
			}
		}

		int t = Ext_ptr[neighbor_index--];
		for (int ix = nx - 1; ix >= 0; --ix)
		{
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(ix, 0, iz);
				if (ix < nx - 1)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
						LocExt_colind[j] = t + nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
					LocExt_colind[j] = t + nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(0, 0, jz - 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
						LocExt_values[j * block_nnz + _B(0, 0, jz) * block_size + _B(bx - 1, by - 1, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(0, 0, jz + 1) * block_size + _B(bx - 1, by - 1, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
						LocExt_colind[j] = t + nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
					}
				}
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(jx, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jx > 0)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(jx - 1, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
							LocExt_values[j * block_nnz + _B(jx - 1, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(jx - 1, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
						}
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(jx, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
						LocExt_values[j * block_nnz + _B(jx, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(jx, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
						if (jx < bx - 1)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(jx + 1, 0, jz - 1) * block_size + _B(jx, by - 1, jz)] = c;
							LocExt_values[j * block_nnz + _B(jx + 1, 0, jz) * block_size + _B(jx, by - 1, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(jx + 1, 0, jz + 1) * block_size + _B(jx, by - 1, jz)] = c;
						}
					}
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
						LocExt_values[j * block_nnz + _B(jx, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
					}
				}
				if (ix > 0)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
						LocExt_colind[j] = t - nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
					LocExt_colind[j] = t - nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, 0, jz - 1) * block_size + _B(0, by - 1, jz)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, 0, jz) * block_size + _B(0, by - 1, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, 0, jz + 1) * block_size + _B(0, by - 1, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(ix, ny - 1, iz)];
						LocExt_colind[j] = t - nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
					}
				}
			}
		}

		if (pz > 0)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int ix = nx - 1; ix >= 0; --ix)
			{
				Ext_map[--t] = _N(ix, 0, nz - 1);
				if (ix < nx - 1)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, 0)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
				}
				int j = --LocExt_rowptr[_N(ix, ny - 1, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					if (jx > 0)
						LocExt_values[j * block_nnz + _B(jx - 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
						LocExt_values[j * block_nnz + _B(jx, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
					if (jx < bx - 1)
						LocExt_values[j * block_nnz + _B(jx + 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
				}
				if (ix > 0)
				{
					int j = --LocExt_rowptr[_N(ix, ny - 1, 0)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
				}
			}
		}
	}

	if (pz < Pz - 1)
	{
		int t = Ext_ptr[neighbor_index--];
		for (int ix = nx - 1; ix >= 0; --ix)
		{
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(ix, iy, 0);
				if (ix < nx - 1)
				{
					if (iy < ny - 1)
					{
						int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
						LocExt_colind[j] = t + ny + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, 0) * block_size + _B(bx - 1, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
					LocExt_colind[j] = t + ny;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(0, jy - 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(0, jy, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(0, jy + 1, 0) * block_size + _B(bx - 1, jy, bz - 1)] = c;
					}
					if (iy > 0)
					{
						int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
						LocExt_colind[j] = t + ny - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
					}
				}
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
							LocExt_values[j * block_nnz + _B(jx, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, 0, 0) * block_size + _B(jx, by - 1, bz - 1)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					for (int jy = 0; jy < by; ++jy)
					{
						if (jx > 0)
						{
							if (jy > 0)
								LocExt_values[j * block_nnz + _B(jx - 1, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							LocExt_values[j * block_nnz + _B(jx - 1, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							if (jy < by - 1)
								LocExt_values[j * block_nnz + _B(jx - 1, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
						}
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(jx, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(jx, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(jx, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
						if (jx < bx - 1)
						{
							if (jy > 0)
								LocExt_values[j * block_nnz + _B(jx + 1, jy - 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							LocExt_values[j * block_nnz + _B(jx + 1, jy, 0) * block_size + _B(jx, jy, bz - 1)] = c;
							if (jy < by - 1)
								LocExt_values[j * block_nnz + _B(jx + 1, jy + 1, 0) * block_size + _B(jx, jy, bz - 1)] = c;
						}
					}
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
							LocExt_values[j * block_nnz + _B(jx, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
					}
				}
				if (ix > 0)
				{
					if (iy < ny - 1)
					{
						int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
						LocExt_colind[j] = t - ny + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
					LocExt_colind[j] = t - ny;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, jy, 0) * block_size + _B(0, jy, bz - 1)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
					}
					if (iy > 0)
					{
						int j = --LocExt_rowptr[_N(ix, iy, nz - 1)];
						LocExt_colind[j] = t - ny - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
					}
				}
			}
		}
	}

	if (pz > 0)
	{
		int t = Ext_ptr[neighbor_index--];
		for (int ix = nx - 1; ix >= 0; --ix)
		{
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(ix, iy, nz - 1);
				if (ix < nx - 1)
				{
					if (iy < ny - 1)
					{
						int j = --LocExt_rowptr[_N(ix, iy, 0)];
						LocExt_colind[j] = t + ny + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, 0, bz - 1) * block_size + _B(bx - 1, by - 1, 0)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, iy, 0)];
					LocExt_colind[j] = t + ny;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(0, jy - 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
						LocExt_values[j * block_nnz + _B(0, jy, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(0, jy + 1, bz - 1) * block_size + _B(bx - 1, jy, 0)] = c;
					}
					if (iy > 0)
					{
						int j = --LocExt_rowptr[_N(ix, iy, 0)];
						LocExt_colind[j] = t + ny - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
					}
				}
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(ix, iy, 0)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
							LocExt_values[j * block_nnz + _B(jx, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, 0, bz - 1) * block_size + _B(jx, by - 1, 0)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(ix, iy, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					for (int jy = 0; jy < by; ++jy)
					{
						if (jx > 0)
						{
							if (jy > 0)
								LocExt_values[j * block_nnz + _B(jx - 1, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							LocExt_values[j * block_nnz + _B(jx - 1, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							if (jy < by - 1)
								LocExt_values[j * block_nnz + _B(jx - 1, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
						}
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(jx, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
						LocExt_values[j * block_nnz + _B(jx, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(jx, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
						if (jx < bx - 1)
						{
							if (jy > 0)
								LocExt_values[j * block_nnz + _B(jx + 1, jy - 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							LocExt_values[j * block_nnz + _B(jx + 1, jy, bz - 1) * block_size + _B(jx, jy, 0)] = c;
							if (jy < by - 1)
								LocExt_values[j * block_nnz + _B(jx + 1, jy + 1, bz - 1) * block_size + _B(jx, jy, 0)] = c;
						}
					}
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(ix, iy, 0)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
							LocExt_values[j * block_nnz + _B(jx, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
					}
				}
				if (ix > 0)
				{
					if (iy < ny - 1)
					{
						int j = --LocExt_rowptr[_N(ix, iy, 0)];
						LocExt_colind[j] = t - ny + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, iy, 0)];
					LocExt_colind[j] = t - ny;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, jy, bz - 1) * block_size + _B(0, jy, 0)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
					}
					if (iy > 0)
					{
						int j = --LocExt_rowptr[_N(ix, iy, 0)];
						LocExt_colind[j] = t - ny - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
					}
				}
			}
		}
	}

	if (py > 0)
	{
		if (pz < Pz - 1)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int ix = nx - 1; ix >= 0; --ix)
			{
				Ext_map[--t] = _N(ix, ny - 1, 0);
				if (ix < nx - 1)
				{
					int j = --LocExt_rowptr[_N(ix, 0, nz - 1)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(ix, 0, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					if (jx > 0)
						LocExt_values[j * block_nnz + _B(jx - 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(jx, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
					if (jx < bx - 1)
						LocExt_values[j * block_nnz + _B(jx + 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
				}
				if (ix > 0)
				{
					int j = --LocExt_rowptr[_N(ix, 0, nz - 1)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
				}
			}
		}

		int t = Ext_ptr[neighbor_index--];
		for (int ix = nx - 1; ix >= 0; --ix)
		{
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(ix, ny - 1, iz);
				if (ix < nx - 1)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(ix, 0, iz)];
						LocExt_colind[j] = t + nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, 0) * block_size + _B(bx - 1, 0, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, 0, iz)];
					LocExt_colind[j] = t + nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(0, by - 1, jz - 1) * block_size + _B(bx - 1, 0, jz)] = c;
						LocExt_values[j * block_nnz + _B(0, by - 1, jz) * block_size + _B(bx - 1, 0, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(0, by - 1, jz + 1) * block_size + _B(bx - 1, 0, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(ix, 0, iz)];
						LocExt_colind[j] = t + nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
					}
				}
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(ix, 0, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(jx, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, by - 1, 0) * block_size + _B(jx, 0, bz - 1)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(ix, 0, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jx > 0)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(jx - 1, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
							LocExt_values[j * block_nnz + _B(jx - 1, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(jx - 1, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
						}
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(jx, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
						LocExt_values[j * block_nnz + _B(jx, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(jx, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
						if (jx < bx - 1)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(jx + 1, by - 1, jz - 1) * block_size + _B(jx, 0, jz)] = c;
							LocExt_values[j * block_nnz + _B(jx + 1, by - 1, jz) * block_size + _B(jx, 0, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(jx + 1, by - 1, jz + 1) * block_size + _B(jx, 0, jz)] = c;
						}
					}
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(ix, 0, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jx = 0; jx < bx; ++jx)
					{
						if (jx > 0)
							LocExt_values[j * block_nnz + _B(jx - 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
						LocExt_values[j * block_nnz + _B(jx, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
						if (jx < bx - 1)
							LocExt_values[j * block_nnz + _B(jx + 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
					}
				}
				if (ix > 0)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(ix, 0, iz)];
						LocExt_colind[j] = t - nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(ix, 0, iz)];
					LocExt_colind[j] = t - nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz - 1) * block_size + _B(0, 0, jz)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz) * block_size + _B(0, 0, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz + 1) * block_size + _B(0, 0, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(ix, 0, iz)];
						LocExt_colind[j] = t - nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
					}
				}
			}
		}

		if (pz > 0)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int ix = nx - 1; ix >= 0; --ix)
			{
				Ext_map[--t] = _N(ix, ny - 1, nz - 1);
				if (ix < nx - 1)
				{
					int j = --LocExt_rowptr[_N(ix, 0, 0)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(0, by - 1, bz - 1) * block_size + _B(bx - 1, 0, 0)] = c;
				}
				int j = --LocExt_rowptr[_N(ix, 0, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jx = 0; jx < bx; ++jx)
				{
					if (jx > 0)
						LocExt_values[j * block_nnz + _B(jx - 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
						LocExt_values[j * block_nnz + _B(jx, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
					if (jx < bx - 1)
						LocExt_values[j * block_nnz + _B(jx + 1, by - 1, bz - 1) * block_size + _B(jx, 0, 0)] = c;
				}
				if (ix > 0)
				{
					int j = --LocExt_rowptr[_N(ix, 0, 0)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
				}
			}
		}
	}

	if (px > 0)
	{
		if (py < Py - 1)
		{
			if (pz < Pz - 1)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(nx - 1, 0, 0);
				int j = --LocExt_rowptr[_N(0, ny - 1, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
			}

			int t = Ext_ptr[neighbor_index--];
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(nx - 1, 0, iz);
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(0, ny - 1, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(0, ny - 1, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jz = 0; jz < bz; ++jz)
				{
					if (jz > 0)
						LocExt_values[j * block_nnz + _B(bx - 1, 0, jz - 1) * block_size + _B(0, by - 1, jz)] = c;
					LocExt_values[j * block_nnz + _B(bx - 1, 0, jz) * block_size + _B(0, by - 1, jz)] = c;
					if (jz < bz - 1)
						LocExt_values[j * block_nnz + _B(bx - 1, 0, jz + 1) * block_size + _B(0, by - 1, jz)] = c;
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(0, ny - 1, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
				}
			}

			if (pz > 0)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(nx - 1, 0, nz - 1);
				int j = --LocExt_rowptr[_N(0, ny - 1, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
			}
		}

		if (pz < Pz - 1)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(nx - 1, iy, 0);
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(0, iy, nz - 1)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(0, iy, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					if (jy > 0)
						LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
					LocExt_values[j * block_nnz + _B(bx - 1, jy, 0) * block_size + _B(0, jy, bz - 1)] = c;
					if (jy < by - 1)
						LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(0, iy, nz - 1)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
				}
			}
		}

		int t = Ext_ptr[neighbor_index--];
		for (int iy = ny - 1; iy >= 0; --iy)
		{
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(nx - 1, iy, iz);
				if (iy < ny - 1)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(0, iy, iz)];
						LocExt_colind[j] = t + nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, 0) * block_size + _B(0, by - 1, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(0, iy, iz)];
					LocExt_colind[j] = t + nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, 0, jz - 1) * block_size + _B(0, by - 1, jz)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, 0, jz) * block_size + _B(0, by - 1, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, 0, jz + 1) * block_size + _B(0, by - 1, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(0, iy, iz)];
						LocExt_colind[j] = t + nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
					}
				}
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(0, iy, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, jy, 0) * block_size + _B(0, jy, bz - 1)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, 0) * block_size + _B(0, jy, bz - 1)] = c;
					}
				}
				int j = --LocExt_rowptr[_N(0, iy, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jy > 0)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, jz - 1) * block_size + _B(0, jy, jz)] = c;
							LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, jz) * block_size + _B(0, jy, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, jz + 1) * block_size + _B(0, jy, jz)] = c;
						}
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, jy, jz - 1) * block_size + _B(0, jy, jz)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, jy, jz) * block_size + _B(0, jy, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, jy, jz + 1) * block_size + _B(0, jy, jz)] = c;
						if (jy < by - 1)
						{
							if (jz > 0)
								LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, jz - 1) * block_size + _B(0, jy, jz)] = c;
							LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, jz) * block_size + _B(0, jy, jz)] = c;
							if (jz < bz - 1)
								LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, jz + 1) * block_size + _B(0, jy, jz)] = c;
						}
					}
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(0, iy, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jy = 0; jy < by; ++jy)
					{
						if (jy > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, jy, bz - 1) * block_size + _B(0, jy, 0)] = c;
						if (jy < by - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
					}
				}
				if (iy > 0)
				{
					if (iz < nz - 1)
					{
						int j = --LocExt_rowptr[_N(0, iy, iz)];
						LocExt_colind[j] = t - nz + 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
					}
					int j = --LocExt_rowptr[_N(0, iy, iz)];
					LocExt_colind[j] = t - nz;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					for (int jz = 0; jz < bz; ++jz)
					{
						if (jz > 0)
							LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz - 1) * block_size + _B(0, 0, jz)] = c;
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz) * block_size + _B(0, 0, jz)] = c;
						if (jz < bz - 1)
							LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz + 1) * block_size + _B(0, 0, jz)] = c;
					}
					if (iz > 0)
					{
						int j = --LocExt_rowptr[_N(0, iy, iz)];
						LocExt_colind[j] = t - nz - 1;
						BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
					}
				}
			}
		}

		if (pz > 0)
		{
			int t = Ext_ptr[neighbor_index--];
			for (int iy = ny - 1; iy >= 0; --iy)
			{
				Ext_map[--t] = _N(nx - 1, iy, nz - 1);
				if (iy < ny - 1)
				{
					int j = --LocExt_rowptr[_N(0, iy, 0)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, 0, bz - 1) * block_size + _B(0, by - 1, 0)] = c;
				}
				int j = --LocExt_rowptr[_N(0, iy, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jy = 0; jy < by; ++jy)
				{
					if (jy > 0)
						LocExt_values[j * block_nnz + _B(bx - 1, jy - 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
					LocExt_values[j * block_nnz + _B(bx - 1, jy, bz - 1) * block_size + _B(0, jy, 0)] = c;
					if (jy < by - 1)
						LocExt_values[j * block_nnz + _B(bx - 1, jy + 1, bz - 1) * block_size + _B(0, jy, 0)] = c;
				}
				if (iy > 0)
				{
					int j = --LocExt_rowptr[_N(0, iy, 0)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
				}
			}
		}

		if (py > 0)
		{
			if (pz < Pz - 1)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(nx - 1, ny - 1, 0);
				int j = --LocExt_rowptr[_N(0, 0, nz - 1)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
			}

			int t = Ext_ptr[neighbor_index--];
			for (int iz = nz - 1; iz >= 0; --iz)
			{
				Ext_map[--t] = _N(nx - 1, ny - 1, iz);
				if (iz < nz - 1)
				{
					int j = --LocExt_rowptr[_N(0, 0, iz)];
					LocExt_colind[j] = t + 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, 0) * block_size + _B(0, 0, bz - 1)] = c;
				}
				int j = --LocExt_rowptr[_N(0, 0, iz)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				for (int jz = 0; jz < bz; ++jz)
				{
					if (jz > 0)
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz - 1) * block_size + _B(0, 0, jz)] = c;
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz) * block_size + _B(0, 0, jz)] = c;
					if (jz < bz - 1)
						LocExt_values[j * block_nnz + _B(bx - 1, by - 1, jz + 1) * block_size + _B(0, 0, jz)] = c;
				}
				if (iz > 0)
				{
					int j = --LocExt_rowptr[_N(0, 0, iz)];
					LocExt_colind[j] = t - 1;
					BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
					LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
				}
			}

			if (pz > 0)
			{
				int t = Ext_ptr[neighbor_index--];
				Ext_map[--t] = _N(nx - 1, ny - 1, nz - 1);
				int j = --LocExt_rowptr[_N(0, 0, 0)];
				LocExt_colind[j] = t;
				BSRBlockFill(block_size, 0.0, LocExt_values + j * block_nnz);
				LocExt_values[j * block_nnz + _B(bx - 1, by - 1, bz - 1) * block_size + _B(0, 0, 0)] = c;
			}
		}
	}

	A.Free();
	A.comm = comm;
	A.LocLoc.size[0] = size;
	A.LocLoc.size[1] = size;
	A.LocLoc.block_size = block_size;
	A.LocLoc.rowptr = LocLoc_rowptr;
	A.LocLoc.colind = LocLoc_colind;
	A.LocLoc.values = LocLoc_values;
	A.LocExt.size[0] = size;
	A.LocExt.size[1] = Ext_ptr[num_neighbors];
	A.LocExt.block_size = block_size;
	A.LocExt.rowptr = LocExt_rowptr;
	A.LocExt.colind = LocExt_colind;
	A.LocExt.values = LocExt_values;
	A.num_neighbors = num_neighbors;
	A.neighbor_ranks = neighbor_ranks;
	A.Ext_ptr = Ext_ptr;
	A.Ext_map = Ext_map;
}

}